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We develop a theoretical approach to study the transient dynamics and the time-dependent statis¬ 
tics for the Anderson-Holstein model in the regime of strong electron-phonon coupling. For this pur¬ 
pose we adapt a recently introduced diagrammatic approach to the time domain. The generating 
function for the time-dependent charge transfer probabilities is evaluated numerically by discretizing 
the Keldysh contour. The method allows us to analyze the system evolution to the steady state 
after a sudden connection of the dot to the leads, starting from different initial conditions. Simple 
analytical results are obtained in the regime of very short times. We study in particular the apparent 
bistable behavior occurring for strong electron-phonon coupling, small bias voltages and a detuned 
dot level. The results obtained are in remarkable good agreement with numerically exact results 
obtained by Quantum Monte Carlo methods. We analyze the waiting time distribution and charge 
transfer probabilities, showing that only a single electron transfer is responsible for the rich struc¬ 
ture found in the short times regime. A universal scaling (independent of the model parameters) is 
found for the relative amplitude of the higher order current cumulants in the short times regime, 
starting from an initially empty dot. We finally analyze the convergence to the steady state of the 
differential conductance and of the differential Fano factor at the inelastic threshold, which exhibits 
a peculiar oscillatory behavior. 


I. INTRODUCTION 

The study of time-dependent current fluctuations in 
nanoscale conductors is of great importance as it can 
provide information on the interactions and quantum 
correlations between electrons While these stud¬ 

ies have been traditionally restricted to the stationary 
regime (corresponding to long measuring times), the ad¬ 
vent of single electron sources has triggered the in¬ 
terest in the short-times behavior. On the one hand this 
knowledge would be useful to fully characterize the single 
electron emitters in the high frequency range @ . On the 
other hand, understanding the short time dynamics is a 
necessary requirement for the use of nanodevices in the 
detection of individual electrons 0 ■ 

This context suggests the need of developing new 
methods to characterize the statistics in the time-domain. 
The concept of Waiting Time Distributions (WTD) is 
well known in the field of quantum optics and stochastic 
processes 0 but has been more recently introduced in 
electronic transport [ 1 ]. Here it has been studied in the 
incoherent r^ime using master equations both within 
Markovian (oMllj and non-Markovian [T^ approxima¬ 
tions. The extension of these studies to the coherent 
(but non-interacting) regime is a quite recent develop¬ 
ment. For this case a scattering approach has been in¬ 
troduced [l^ and adapted later to tight-binding models 
M- Other approach to the problem is provided by non¬ 
equilibrium Green functions methods, which have been 
discussed in Ref. [l^ and applied to analyze the tran¬ 
sient dynamics of non-interacting quantum dots in Refs. 


mill. 

Green function methods are in principle the most 
appropriate to study the effect of interactions in the 
time-dependent statistical properties of quantum coher¬ 
ent conductors. However, their application to this case 
has remained essentially unexplored. In the present work 
we provide some initial steps in this direction by analyz¬ 
ing the Anderson-Holstein model in the polaronic regime. 
This simple model provides the basis to understand quite 
complex non-equilibrium phenomena occurriM in actual 
systems such as phonon-assisted tunneling [I^j and 
Franck-Condon blockade [2^ . 

While the stationary transport properties of the 
Anderson-Holstein model have been extensively analyzed 
[2ll - l^ . their analysis in the time domain has been 
much less analyzed [l^. Recent calculations, based 
on numerically exact methods like diagrammatic quan¬ 
tum Monte Carlo (diagMC) [2^ have indicated that for 
strong electron-phonon coupling there exists a regime in 
which different initial conditions lead to different trans¬ 
port properties at short times. In a subsequent work [2^ 
it was demonstrated that this apparent bistability actu¬ 
ally corresponds to a long transient dynamics leading to 
blocking-deblocking events associated to the polaron dy¬ 
namics. More recent work has confirmed the analy¬ 
sis and also explored the effect of including the dot-leads 
Coulomb repulsion. The approach to the steady state 
has also been analyzed for the Anderson-Holstein model 
including a continuous distribution of phonon modes . 

In spite of these efforts, none of these works have an¬ 
alyzed the time evolution of the noise properties of the 
Anderson-Holstein model as the system approaches the 


2 


steady state. This deficit connects with the above men¬ 
tioned lack of studies of time-dependent statistics for 
interacting systems in the quantum transport regime. 
Unfortunately, numerically exact methods like quan¬ 
tum Monte Car lo 1321 jSSj or Numerical Renormalization 
Group (NRG) have not yet been adapted to noise 
studies. 

To circumvent these difficulties, the present work in¬ 
troduces a generalization to the time-dependent case of 
a simple analytical approach called Dressed Tunneling 
Approximation (DTA), which was shown to give a good 
description of the spectral and transport properties in the 
stationary limit for the polaronic regime |35| . As a first 
step we check that the method provides results for the 
time-evolution of the mean current and the dot charge 
which are in good agreement with numerically exact re¬ 
sults of Ref. . We then study the transient WTD and 
the evolution of the current cumulants, showing that in¬ 
teractions tend to increase the characteristic times for 
relaxation towards the steady state and also enhance the 
asymmetry in the charge transfer probability distribu¬ 
tions. We also study the scaling of the relative ampli¬ 
tudes of the transient cumulants of higher order and find 
a very robust universal behavior, that we demonstrate 
using analytical arguments. Finally, we analyze the con¬ 
vergence to the steady state of the differential conduc¬ 
tance and of the differential Fano factors at the inelastic 
threshold V = wq and also at U = 2wo. 

II. MODEL AND BASIC THEORETICAL 
FORMULATION 

We consider the simplest spinless Anderson-Holstein 
model in which a single electronic level is coupled to a lo¬ 
calized vibrational mode. Electrons can tunnel from this 
resonant level into a left (L) and a right electrode (R). 
We shall generically refer to this central region, which 
can represent either a molecule, an atomic chain or a 
quantum dot, as the “dot” region. The corresponding 
Hamiltonian is given hy H = Hi^ads + H^ot + with 
(in natural units, h = ks = & = rrie = 1) 

Hdot = [eo + A (a'l' -I- a)] d^d + luq a^a , (1) 

where eo is the bare electronic level, A is the electron- 
phonon coupling constant and luq is the frequency of the 
localized vibration. The electron (phonon) creation op¬ 
erator in the dot is denoted by d^ (a'l'). On the other 

hand, Hieads = Y.vk ^vkclf,Cvk corresponds to the non¬ 
interacting leads Hamiltonian {u = L, R) where e^k are 
the leads electron energies and are the corresponding 
creation operators. The bias voltage applied to the junc¬ 
tion is imposed by shifting symmetrically the chemical 
potential of the electrodes V — — k'R- 

The tunneling processes are described by 


where 'ji.k are the tunneling amplitudes. 

To address the polaronic regime we perform the so- 
called Lang-Firsov unitary transformation which 
eliminates the linear term in the electron-phonon cou¬ 
pling [13 

= = g = —. (3) 

Wo 

Using this transformation 

Hdot = ed^ d + LUoa^a , (4) 

where e = cq — A^/wq. The tunneling Hamiltonian is 
transformed as 

Ht = '^ (j,.k clf, Xd + h.c.) , (5) 

uk 

where X = exp [^(a — a^)] is the phonon cloud operator. 
On the other hand, the free leads Hamiltonian remains 
invariant. For later use it is useful to introduce the tun¬ 
neling rates Fi, = Im^^ |7i/fcP/(w —jO’*' —Ci/fc) which are 
approximated by constants in the so-called wide band 
approximation, and denote F = F/, -|- F/j. 

In the present work we focus on the transient dynamics 
which corresponds to the evolution of the system from an 
initial t = 0 state when the dot is suddenly connected to 
both leads. The corresponding statistical properties of 
the transferred charges can be obtained from the gener¬ 
ating function (GF) 

OO 

Z{x.t)= Y. ( 6 ) 

q— — co 

where Pq{t) denotes the probability of transferring q 
charges through the dot in the measuring time t. The 
GF is in turn related to the Cumulant Generating 
Function (CGF) by <S'(x,t) = logZ(x,t), which gen¬ 
erates the time-dependent charge cumulants Ck(t) = 
{—i)'^d^S/dx'^\x=o- We further define ((/^(t))) = 

dCk {t)/dt, which tend to the zero-frequency steady state 
current cumulants when t —>■ oo. Another quantity of in¬ 
terest to characterize the transient statistics is the wait¬ 
ing time distribution W{t). This can be related to the 
so-called idle time probability n(t) = Po{t), defined as 

H 

= Az(x,i). (7) 

While in a stationary situation the definition of the 
WTD requires a two-time measurement [l3l| , in the tran¬ 
sient case the initial time is fixed at t = 0 and one can 
define a single time measurement WTD as 1,13 


Ht = Y ^Ik d > ( 2 ) 

uk 


w{t) 


dll{t) 

dt 


( 8 ) 
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which gives the probability that the first electron is de¬ 
tected at a certain time t after the connection to the 
leads. 

The GF can be written [T^ as an average of the evo¬ 
lution operator over the Keldysh contour, shown in Fig. 

ID 


Z{x,t) =< Tcexp{-i J > , (9) 

where is the system Hamiltonian with a counting field 
x{t) which take the values ±x on the two branches of the 
Keldysh contour entering as a phase factor modulating 
the tunnel Hamiltonian, i.e. 

= E 4 Xd + h.c.) . (10) 

yk 

Notice that different charge and current cumulants can 
be defined depending on how the phase x(t) is distributed 
on the left and on the right tunnel couplings. For in¬ 
stance, taking xl = x(t) and XR — Oj generates 

the current and charge transfer cumulants through the 
interface between the left lead and the dot. This is the 
choice that we shall select for the rest of the paper, unless 
specified explicitly. 


A. Non-interacting case 

In Refs.[l^, it has been shown by path-integral 
methods that in the non-interacting case Z(x,t) can be 
expressed as the following Fredholm determinant, defined 
on the Keldysh contour 


Z(x,t) =det 




( 11 ) 


where G and G denote the dot Keldysh Green functions, 
go corresponds to the uncoupled dot case and E are the 
self-energies due to the coupling to the leads. In the 
quantities G and E the tilde indicates the inclusion of 
the counting field in the tunnel amplitudes. 


t=0 




FIG. 1. Keldysh contour considered to analyze the transient 
regime, x indicates the counting field changing sign on the 
two branches of the contour and At corresponds to the time 
step in the discretized calculation of the generating function 

z{x,t)- 


As shown in a simple discretized version of the 
inverse free dot Green function on the Keldysh contour 
is 


(t 

h_ -1 

-p\ 


I 

-1 



h_|_ —I 


V 

h+ -1/ 

2Nx2N 


( 12 ) 

where h± = I ieoAt, At indicates the time step in 
the discretization with N = t/At. In this expression p 
determines the initial dot charge nj, by Ud = p/{l + p). 

On the other hand, the self-energies are given by 


E“^(t,t') =a/30(l)0(i')E4e*(“-^)^^/"4(<,t') , 

yk 

(13) 

where g‘^^{t,t') = -i{Tcc^k{to,)cl^{t'p)), with a,/3 = 
-h, — are the Keldysh Green functions of the uncoupled 
leads. In Appendix lAl we discuss the discretization pro¬ 
cedure and give the explicit expressions for these self¬ 
energies in the discretized contour. 

It should be noticed that while E are 2x2 block 
Toeplitz matrices depending only on the time arguments 
difference, g^^ deviates from a perfect block Toeplitz ma¬ 
trix due to the [N -\- 1, N) and (I, 2N) entries associated 
with the closing of the Keldysh contour and the initial 
condition respectively. The connection with the theory 
of Toeplitz determinants is an interesting issue that 
goes beyond the scope of the present work and will be 
discussed elsewhere. 


B. Interacting case: Dressed tunneling 
approximation 


We now discuss the generalization of the theory to the 
interac ting case within the approximation introduced in 
in Ref. [3^. For this purpose we start from the counting 
field functional derivative of the GF 
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p = [ dhY, iLk < Tc - xHh)e-^^^*^'>/^CLkiti)d{ti)) > . (14) 

X Jc f. 

This can be related to the three point Green function (TcX(t)c^^(ti)(i(t')), whose equation of motion is 

{^dt, - eLk) {TcX{t)cl^{t^)d{t')) = , (15) 

which can be integrated yielding 

{TcX{t)cl^{h)d{t'))=-iLk f dt2e-*^(‘^)/"5Lfc(l2,li)(TcX(<)Xt(<2)d^(ti)d(t')) • (16) 


Within the decoupling procedure corresponding 
to the DTA one has {TcX{t)X^t 2 )d^ {ti)d{t')) ~ 

{TcX{t)X^t 2 )){Tcd''^{ti)d{t')) which finally allows us to 
write 


^ = -y dh j dt2Trif |^^^(ti,t2)G(t2,li)| , 

(17) 

where Tr^ denotes trace over the 2x2 Keldysh space 
and f^DTA is the DTA self-energy whose components are 
given by 


energies in the non-interacting case given by Eq. m- 
On the other hand, the propagator A°‘^(1:, t') will be eval¬ 
uated assuming equilibrated phonons (see Appendix IbI) . 

Integrating Eq. (HZl) and imposing the condition 
Z(0, t) = 1 one arrives to the same expression for Z(x, t) 
as in Eq. (EU but replacing S and G by Tjdta and 
Gdta- All the effect of interactions are thus encoded 
in the DTA self-energy Y^dta- More details on the ap¬ 
proximation are given in Appendices [B] and [C] It should 
be noted that this simple structure is valid within DTA 
but in a more general approximation vertex corrections 
would prevent the counting field integration leading to 
Eq. (HII). 


SSTA(i.O = S“''(t,<')A“^(t,t'), (18) 

C. Tunnel and short time limits 

with A“^(t, t') = (TcA(t)Xl(t')) being the phonon cloud 

propagator. In this expression S“^ denote the self- To the lowest order in T we have the expansion 

I 


■^(X)l) —1 + J dti J dt2Trx I ^Sdta(1i, I 2 ) — S£ita(1i, 12 )^ ffo(l2, li)| 


(19) 


which reduces to 

Zix,t) - ^ + dh dt2 (s+OT^(ti,t2)5o“’^(^2,ti) - 1)-f S^+ta(^i>^ 2)5(('“(12,D) (e"*^ - 1)) (20) 

It is interesting to notice that the T —>■ 0 and the t —)• 0 limits should coincide, i.e. the short time behavior is well 
described by the expansion to the lowest order in Ft. This allows us to obtain the short time limit of the GF as 

Z{x,t) ^1 + |(e*^ - 1 )Aloi( 1)[1 - «d] + (e“*^ - l)Aiio(t)nd| , (21) 


with 


ALOlit) = — X 

TT 


E 

n—0 


1 - cos[(a; - e - noJo)t] 
an I duJ -^-r-- JlIU}) 


-w 


(uj — i — nojo)^ 


AlioW = — X 

TT 


E 

n—O 


f ^ ^ - cos[{uj - e + nuJo)t] 
an / duj - - -—--[/L(a;) - 1] 

-w (uj-e + nujo) 


where fLi^j) is the Fermi distribution at the left electrode 
and, at zero temperature, a„ = e“® g'^^jnl. The physical 
interpretation of the amplitudes A^oi and A^iq is trans¬ 
parent in the t —>■ 00 limit where A^Qi/t and A^Qi/t tend 
to the Fermi Golden rule rates derived for the sequential 
tunneling regime (see e.g. Ref. EH). 


{uj — e + nwo)^ 
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As expected in the short time limit Ft 1, the GF in 
Eq. (I?n) involves charge transfer of a single electron, with 
only Po{t), P-i{t) and Pi{t) having a significant weight, 
and are given by 


Pq{t) = ALoi(t)[l — nci]Sq^i — ALio{t)ndSq-i , (22) 

I 


and Po{t) = 1 — Pi{t) — P-i{t). Correspondingly, the 
WTD is proportional to the left current I Lit), he. 

Wit) = -j^Po{t) = \^-^\. (23) 

At zero temperature and neglecting contributions from 
the band edges (which is justified in the wide band ap¬ 
proximation) we obtain 


II = rL(l - 2nd) H- - ctn {Si [{^il + nujo - e)t] (1 - Ud) - Si [{(Ml - nuJo - e)t] Ud} (24) 

TT ^' 

n 


where Si denotes the sine integral function. 

Due to the property lima;_>o Si(a;) = 0 we find that the 
initial current is II = iF/,, with the sign depending on 
the initial charge (ud = 0 or = 1). This property fixes 
also the initial value of the WTD as Wit) = abs (/l( 0)). 
It should be noticed that we are neglecting the system 
evolution on time scales smaller than the inverse of the 
leads bandwidth (see Appendix , which explains why 
the initial current can be non-zero. However, for the 
symmetrized current 1)1) = iVnlhit) — Fi//i(t))/F the 
initial value is zero and the first charge transfer cumulant 
is a continuous function starting from zero regardless of 
the initial condition. 


III. RESULTS 

A. Evolution of mean current and charge: 
comparison to diagMC results 

We start by analyzing the transient behavior of the 
Holstein model for different initial conditions. In Ref. 

it was shown that for certain parameters range the 
model exhibits an apparent bistable behavior. Further 
analysis [H, [s^ revealed that this apparent bistability 
was caused by a long transient regime associated with 
the slow polaronic dynamics. The comparison of the 
symmetrized current and the charge evolution obtained 
within DTA with the numerically exact results obtained 
by diagMC are shown in Figs. [2] and [S] These results 
corresponds to the case of a deep level e = —IOF and 
large phonon frequency ujq = 8F where the apparent 
bistable behavior is more pronounced. There is a re¬ 
markable agreement between the DTA and the diagMC 
results for the current in the U = 5F case. The agreement 
with the current is somewhat poorer for V = 25F, which 
reflects a limitation of the DTA for describing the spec¬ 
tral density around the Fermi energy in this regime (see 
[ 3 ^) but there is a very good agreement in the evolution 
of the mean charge for this case (see Fig. IS. This repre¬ 
sents an improvement with respect to the approximation 


used in Ref. [^ . 



^ ^ - -- - 

‘0.0 0.5 1.0 , 1.5 2.0 2.5 

tr'] 


FIG. 2. (Color online) Symmetrized transient current I(t): 
comparison between DTA (full lines) and diagMC (dots) for 
initially empty (red) and initially full (blue) dot. e = — lOF, 
g = 2, Wo = Sr and F = SF (upper panel) and V = 26r 
(lower panel). The inset illustrates the convergence to the 
steady state for the case V = SF. 


The apparent bistable behavior at short times corre¬ 
sponds actually to a long transient dynamics, as illus¬ 
trated by the DTA results on a longer time scale (see 
inset in Fig. HD. 
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FIG. 3. (Color online) Time-dependent dot occupation: com¬ 
parison between DTA and diagMC for initially empty and 
initially full dot. Same parameters as in Fig. (lower panel). 


B. Waiting time distribution and transient 
statistics 

Further insight on the transient dynamics is provided 
by analyzing the evolution of the WTD and the higher 
current cumulants, see Fig. S) 



FIG. 4. (Color online) Upper panel: waiting time distribution 
for the initially empty (red) and initially occupied (blue) cases 
for the same parameters as in Fig. [2] (lower panel). Inset: 
probabilities Po (full line) and Pi (dashed line) for the initially 
empty state. The lower panel shows the transient dynamics 
of current cumulants with k = 1,2, 3, 4, 5,6 (from 

top to bottom at short times) for the initially empty case. 

As can be observed, in the initially empty case the 
WTD exhibits a non-monotonous decrease with small 
steps at time ~ n2'K jujQ associated with the polaron dy¬ 
namics. The current cumulants in this short time regime 
have an increasing amplitude with increasing cumulant 
order (we come back to this point below). The inset 
in the upper panel of Fig. 0] indicates that this short 
time dynamics is associated with a single electron trans¬ 
fer, with only Pq and Pi being non-negligible. Thus the 
waiting time distribution follows essentially the current 
at short times. On the other hand, for the initially occu¬ 


pied state the current dynamics is almost blocked in this 
short time regime. 

All the previous results indicate that the typical times 
for the transient regime are increased by the effect of in¬ 
teractions. A more clear picture of this effect is provided 
by Fig. [5] where W{t) is shown for increasing values of g 
for the initially empty and initially occupied states. 

At short time scales {t < F”^) and for weak cou¬ 
pling to the leads, in the initially empty case, the WTD 
is well approximated by Eq. (IMl) . Thus, it exhibits 
an initial linear behavior fixed hy W{t) w rL/2 -|- 
2 Fi/7r —e — nuio)t followed by an extremum at 

t ~ '^{\J2n<^nigL-e-nUJo)/J2n “ U-Wq)^ | 

In the non-interacting case g = 0 only the n = 0 term 
contributes and the WTD exhibits an initial positive 
slope for the parameters in Fig. [5l reaching a maximum 
peak located at t ~ 7r/|U/2 — e|. Increasing the coupling 
strength, the initial slope decreases and becomes even¬ 
tually negative. For g « 1.3, the peak becomes a dip 
indicating the transition into the strong coupling regime. 
The probabilities Pi(t) and P 2 it) shown in the insets al¬ 
low to visualize the injection of the first and second elec¬ 
tron, and their ralentization with increasing interaction. 



FIG. 5. (Color online) Waiting time distribution for increas¬ 
ing g (0 (green), 1 (blue) and 1.5 (red)) for initially empty 
(upper panel) and initially full (lower panel) with e = —F, 
ujo = 2r and U = 5F. The insets show the corresponding 
probabilities Pi{t) and P 2 {t) (upper panel) and Pi(t) and 
P-i{t) (lower panel). 

On the other hand, the initially occupied case ex¬ 
hibits a different short time linear scaling for the WTD 
W{t) ss ri/2 - ‘2rL/7rJ2jiOiniV/2 - e + nwo)^ followed 
by a dip at very short times t rs 7r{| — e + 

^n{gL — e + which is associated to 

the blocking effect of the occupied level. Contrary to the 
empty case, the evolution of the dip is monotonous with 
the interaction strength: the dip depth increases and its 
position shifts towards smaller times with increasing g. 
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As shown by the insets, the backward flow probability 
P-i{t) is quite significant in this case, contrary to the 
initially empty case. 

Similarly, the increase of g has an impact in the evo¬ 
lution and the stationary limit of the higher order cumu- 
lants ((7^)). These are shown in Fig. [6]for fc = 2, 3 and 
4. The cumulants are normalized to the time dependent 
mean current, which allows to appreciate more clearly 
the differences with increasing g. As a general feature, 
both the interacting and non-interacting cases exhibit an 
increase in the transient amplitude with increasing cumu- 
lant order, as already mentioned in connection to Fig. U) 
The effect of increasing g is twofold: first, it slows down 
the dynamics and second, the relative asymptotic values 
of the cumulants are larger than the non-interacting ones. 
In the initially occupied case (lower panels in Figl6]) the 
same effects can be observed. The divergent relative am¬ 
plitudes at very short times are due to the change of sign 
of the mean current occurring in this case. 



FIG. 6. (Color online) Higher order cumulants (normalized to 
the mean current) for increasing g values (same color code as 
in Fig. ID) for the initially empty (upper panels) and initially 
full cases (lower panels). Notice that in the lower panels the 
time axis starts at times slightly larger than 0 in order to 
avoid the divergent behavior of the normalized cumulants in 
the f ^ 0 limit. 


C. Universal scaling of normalized transient 
cumulants 

In spite of the renormalization of the characteristic 
times and of the asymptotic cumulants introduced by 
interactions, for the initially empty case some “univer¬ 
sal” features can be identified. For instance, as shown in 
the lower panel of Fig. [3 the maximum amplitude of the 
normalized transient cumulants Ck(t)/Ci{t) are found to 
follow a scaling which is independent of the value of the 
interaction parameter g. This scaling is also found to be 
extremely robust with other model parameters like e or 
Wq. 

We can explain this behavior by considering the short 
time limit discussed in Sect. Ill Cl For the initially empty 
case gQ~{t,t') — 0 and as Eq. (1^ indicates, the current 
flow is unidirectional. Moreover, in this limit all deriva¬ 
tives of the GF are equal, i.e. (—f)”9"Z(t, x)/9x"Jo = 
x{t), corresponding to the situation where just a single 
electron is involved, as shown in the inset of Fig. 21 Cor¬ 
respondingly, the CGF at short times can be written as 

S'(x,t) ~ log [l-ba;(t) - I)] . (25) 

A simple Taylor expansion allows to identify the charge 
cumulants as 


= (26) 
9=1 


where u{t) = x{t)/{1 — x{t)). Taking the continuous limit 
and evaluating the above summation as an integral (valid 
at sufficiently large k) we find 


Ck = 


2{k — I)! cos 


k arctan 


(isin)] 


(log^ u + 


kjl 


(27) 


where cr„ = sign(log u). 

While the cosine factor in this expression is respon¬ 
sible for the oscillatory behavior of the cumulants at 
short times, their maximum amplitude is controlled by 
the (k — I)! factor and the power law in the denomi¬ 
nator. As in the region of maximum amplitude typi¬ 
cally log^ u <C TT^, it is straightforward to show from Eq. 

that the cumulants maximum amplitude scale as 
{k — l)!7r“^. This law describes with accuracy the scal¬ 
ing of the relative current cumulants already shown in 
the lower inset of Fig. 0 

The factorial increase of the transient cumulants has 
been already pointed out and demonstrated experimen¬ 
tally in Ref. [^. In contrast to the present study, they 
considered the sequential tunneling regime. However, as 
we show in the upper panels of Fig. [3 the short time be¬ 
havior of the transient cumulants predicted by Eq. (123 
remarkably agrees with the results of Ref. [d^. The 
reason for this agreement is the universality of the gener¬ 
ating function in the short time regime for unidirectional 
transport corresponding to the initially empty dot case. 


































At longer times the results from Ref. [l^ deviate from 
the predicted behavior by Eq. (HZD by a global extra 
exponential decay, which can be associated with the dif¬ 
ferences with the setup considered in Ref. [4l |. 



FIG. 7. (Color online) Upper panels: normalized transient 
current cumulants obtained from Eq. m (fu ll lines) and the 
corresponding results adapted from Ref. [4^ (dashed lines) 
for fc = 4, 5, 6, 7 (left panel) and k = 8,9, 10,11 (right panel). 
For the comparison we have assumed that Ci{t) follows a 
simple law A(l —e“'"*), which sets the time scale. Lower panel: 
scaling of the maximum amplitude of the relative cumulants 
CkjCi (indicated by max(fc) in the figure) as a function of 
k. The relative cumulants are normalized with (fc — 1)! and 
the linear slope in the logarithmic scale indicates scaling as 
{k — l)!7r~*^ (see text). The different symbols correspond to 
the cases ^ = 0, 1 and 1.5 in Fig. [5] 


D. Conductance and Fano factor dynamics at 

V = noJo 

It is also worth analyzing the current and noise dynam¬ 
ics for bias voltages close to the conditions V nojQ. The 
behavior of the stationary noise at the inelastic threshold 
V = LOo has been analyzed in several works in the limit of 
weak interaction [431 - 14^ . showing that it can either ex¬ 
hibit an increase or a decrease due to the opening of the 
inelastic channel. In contrast, for V = 2ujo, the analysis 
of stationary noise in the polaronic regime presented in 
Ref. [ 3 ^ indicates that it exhibits a suppression associ¬ 
ated to the opening of a side band (new elastic channel). 
The transient conductance for V ^ loq and V ^ 2uJo 
are shown in Figs. [8] and [9] for e = 0, g = 1.5 and 
different values of T. For comparison we show the pre¬ 
diction for the tunnel limit, given by Eq. (m as dashed 
lines. To illustrate the behavior of the noise we choose to 


represent the differential Fano factor dF{t)/dV, where 
F{t) = {{P))/{2{{I))), which allows to see more clearly 
the differences in the T —>■ 0 limit. 



FIG. 8. (Color online) Transient conductance (upper panel) 
and differential Fano factor (lower panel) for V = ljq, g = 1.5, 
e = 0 and different values of F/wo (0.05 (blue), 0.25 (green), 
1 (yellow) and 2 (red)). The dashed line in the upper panel 
corresponds to the tunnel limit analytical result of Eq. (PI|). 

As can be observed in these plots, the behavior of the 
transient conductance for T —>■ 0 is remarkably different 
in the two cases: while for V = ujq the conductance ex¬ 
hibits a sequence of up and down steps at t ~ 2n7r/a;o; 
for V = 2wo the sequence corresponds to steps up only. 
The approximate behavior is well captured by the ana¬ 
lytical expression of Eq. (El shown as dashed lines in 
Figs. [5] and H When T increases the step structure in 
the conductance is progressively damped. On the other 
hand, the noise exhibits an interesting different evolu¬ 
tion with increasing T. While for T —>■ 0 the differential 
noise and the conductance are approximately equal (as 
expected for the tunnel limit) for larger T the differen¬ 
tial Fano factor converge to either positive or negative 
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values for V = ujq, or systematically to negative values 
for V = 2ujq. This is consistent with the predictions of 
Ref. [ 3 ^ for the stationary case. It is interesting to re¬ 
mark the difference between the present calculations and 
those of Ref. [4l|. In that work a sequential tunneling 
approach was used, including the effect of the induced 
nonequilibrium phonon population leading to giant Fano 
factors associated to the Frank-Condon blockade effect. 


the simplest case of a flat density of states (wide band 
approximation) within an energy range in the interval 
{—W,W). The leads bandwidth W is assumed to be 
much larger than the tunneling rates Fi,. The discretized 
time-dependent self-energies correspond to the Fourier 
transform of these energy dependent self-energies evalu¬ 
ated at the discrete time mesh defined by tj = which 
at zero-temperature yields 


IV. CONCLUSIONS AND OUTLOOK 

In this work we have presented an analysis of the time- 
dependent statistics of electron transport through a res¬ 
onant level coupled to a localized vibrational mode. We 
have restricted our analysis to the transient regime which 
is established after a sudden connection of the level to the 
leads. For this analysis we have adapted the recently de¬ 
veloped DTA decoupling scheme [35l |. which provides a 
good description of the stationary transport properties in 
the polaronic regime. In spite of its approximate charac¬ 
ter, our analysis has revealed several features of general 
validity. In the first place it shows that interactions tend 
to increase exponentially the relevant time scales for the 
transient dynamics, leading to an apparent bistability at 
short times in certain parameters regime, in agreement 
with a previous analysis by some of us [29|. Second, we 
have demonstrated that in the short time scale and for 
an initially empty dot the higher order cumulants exhibit 
oscillations with a universal scaling amplitude. This uni¬ 
versal character arises due to the fact that a single elec¬ 
tron transfer controls the transport properties at these 
short time scales. Our analysis has furthermore revealed 
a peculiar oscillatory convergence of the conductance and 
the Fano factor at the inelastic threshold V = ujq. 

The present work constitutes a first step in the study 
of time-dependent statistics for the interacting quantum 
coherent transport regime. We envisage several possible 
extensions of this work, like the study of waiting time 
distributions in the stationary regime for interacting sys¬ 
tems; analyzing the effect of unequilibrated phonons and 
extending the study to superconducting systems. 
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Appendix A: Time discretization procedure 

We describe in this appendix the time discretization 
procedure. Por simplicity we discuss here first the non¬ 
interacting case. For describing the leads we consider 




v=L,R 


%* = - E 


^ ii- j) 


u=L,R 


(* - J) 


(Al) 


where > 1. The other self-energy components 
are given by E++ = -0[i- j] E“+ - 0[j - i] and 
E“.“ = -e[i- j\ E+.“ -0[j - i] E“.+ , where 0 [n] is the 
discrete form of the Heaviside function. These last ex¬ 
pressions are not unambiguously defined for i = j. We 
have found that the more stable algorithm corresponds 
to the choice 0 [0] = 1/2. 

The finite flat bandwidth model at zero temperature 
produces current cumulants increasing linearly from zero, 
followed by an oscillatory pattern on the time scale 1/VF, 
produced by the exponential term in Eq. (EH)- This pat¬ 
tern is absent in the symmetrized current cumulants. In 
order to resolve these features, the time step has to be 
smaller than this typical time, which leads to the con¬ 
dition of At < l/kF for a stable algorithm. In order to 
avoid features associated with the finite bandwidth it is 
convenient to generalize the calculation to finite temper¬ 
atures, performing an expansion of the self-energies in 
Matsubara frequencies. In this case, the limit kF —)• 00 is 
well defined, avoiding divergences at time equal to zero. 
The expressions for the self-energies at finite temperature 
are 

v=L,R 

= 2* E - 5[j - k]) , (A2) 

u=L,R 


where the Fermi function can be computed as 




jk 


i'^Rn 0[j - k]( 

n—0 

-e[k - 


,fin{j-k)At 




^[j - k] 

2 ■ 

(A3) 


Here /3„ and Rn represent the poles and the residues of 
the Matsubara expansion. The convergence speed can be 
improved by using the approximated poles and residues 
proposed by T. Ozaki and computed using a con¬ 
tinued fraction. Differently from the zero-temperature 
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finite bandwidth case, the fast oscillatory term is absent, 
and the maximum of the time step is controlled by the 
r parameter. We have found that At < l/(10r) is 
sufficient to warrant the convergence of the numerical 
algorithm. 


In all the calculations presented in this 
work we have introduced a small temperature of the or¬ 
der of O.ir which helps to stabilize the numerical calcu¬ 
lations but does not produce significant deviations from 
the zero-temperature results. 


In the interacting case, the DTA self-energy, given by 
Eq. m, can be written as an infinite sum over polaronic 
weights (see Appendix |B]) . The number of terms needed 
to be added for convergence is higher when increasing the 
electron-phonon coupling, meaning that higher sidebands 
become more important. Considering that the series can 
be truncated with precision enough at a given integer 
value, n, the condition to converge the calculation can 
be obtained following the same reasoning done in the 
non-interacting case, i.e. the polaronic terms generate 
oscillations of a minimum period l/(na;o), which have 
to be resolved for convergence, leading to the condition 
At < l/(na;o). 


Appendix B: Dressed Tunneling Approximation 

In this appendix we summarize the main expressions 
of the Dressed Tunneling Approximation (DTA) used in 
this work. The DTA approximation is built from the 
Dyson equation, where the leads self-energies have been 
dressed with the phonon cloud [s^ 


G = 5o + 5 oE_dtaG, (B1) 

where go corresponds to the undressed dot, is the 

DTA self-energy (as defined in Eq. (IT51) 1 and the " is used 
to denote Keldysh structure. In this expression integra¬ 
tion over internal times is implicitly assumed. Finally, 
dressing again the full Green function, we find the final 
expression as 

Gf^^{t,t') = G^f^{t,t'Wf^{t,t'), (B2) 

where the phonon cloud propagator A“^(t,f') is evalu¬ 
ated assuming equilibrated phonons, i.e. 


OO 

A+-(t,f') = (A-+(f,t'))* = Y. (B3) 

n— — co 


with 


{2g^^np{l + nS^ (B4) 

In being the modified Bessel function of the first kind, 
which is symmetric in the k argument (In = I-n) and 
Up is the Bose factor 1/ (e^““ — l) with /3 = 1/T. The 
remaining components A++(f,f') and A {t,t') are de¬ 
termined by A++(t,f') = 6{t — t')A ^{t,t') + 9{t' — 
t)A^ {t^t') and A {t^t') = 9{t' — t)A -I- 9{t — 


Appendix C: Single pole approximation 


In this Appendix we discuss an approximation that can 
be used in order to compute the average current and pop¬ 
ulation of the level more efficiently. This approximation 
is useful in the regime where the electron-phonon cou¬ 
pling is strong and the evaluation of the Fredholm deter¬ 
minant m becomes computationally more demanding. 
For the evaluation of the current and dot population the 
counting field is not needed and the transformation to 
the triangular representation in the Keldysh formalism 
can be performed [1^. Then, the Dyson equation for the 
retarded component of the Green function is given simply 
by 

G^{t,t') = g^{t,t') + [ dt2 K^{t,t2)G^{t2,t'), 

Jo 

(Cl) 

where K[t — t 2 ) = f dtigo{t — ti)E(fi — t 2 ) is the Kernel 
of the integral equation. This kernel is time translational 
invariant, meaning that it only depends on the difference 
between time arguments, which implies that the retarded 
Green function, solution to the equation, preserves this 
symmetry. In this special case, Eq. m can be Laplace 
transformed arriving to 


G%) 


9o{s) 


(C2) 


The issue of inverting the Laplace transform of a given 
function is in general not simple, because it usually in¬ 
volves the integration in a Bromwich contour of the func¬ 
tion with a given collection of poles. However, if T <C wq, 
the function has only a single dominant pole, r, and the 
retarded Green function of the dot can be written simply 
as 


G^{t - t') ss 9{t - <')Res(G^, s = r) ^ (C3) 


where Res indicates the residue at s = r. This function 
is exponentially decaying, with a decay rate given by the 
real part of r. 

For the calculation of the current and dot charge evo¬ 
lution we need the G^ component, which we analyze 
below for the initially empty case. For the initially oc¬ 
cupied case, similar expressions can be derived by only 
changing the Green function component to the G ^ one. 
The G'^ component can be computed using the kinetic 
equation 

GoTyl(^>^0= [ dtidt2G^{t,ti)Y.^':^^{ti,t2)G^{t2,t'), 
Jo 

(C4) 
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where the 'Eota is the self-energy in our DTA approxi¬ 
mation. The evolution of the average charge is given by 
the imaginary part of this Green function, i.e. nd{t) = 
Im[G^y^(t, t)]. 

Finally, when the coupling strength to both electrodes 
are equal (Fl = Fi?), the symmetrized current can be 
computed as 


where 


m) = 


[ dti 
Jo 


X [TLfLit - h) - TRfRit - ti)]} , (C5) 


= [A+ 

+ (A+-(t,ti) - A-+(t,ti))G+y^(t,ti)] . 


(C6) 
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